敏感性分析之Stata实操：控制变量内生时的系数敏感性分析-regsensitivity
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2022 暑期班
作者：陈卓然 (中山大学)
邮箱：chenzhr25@mail2.sysu.edu.cn
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 引言
2. 方法初探
3. 边界
4. 截断前沿
5. 回归结果的解读
6. 参考文献
7. 相关推文
1. 引言
因果效应识别的最大挑战之一就是遗漏变量。在线性模型中，当一个控制变量与遗漏变量相关时，这个控制变量的系数将是有偏的。为看清这一点，不妨回忆一下最简单的二元线性回归模型。假定真实模型如下：
从而正确的估计结果为：
然而由于某种原因，我们实际估计的模型如下：
从而错误的估计结果为:
我们将式 (1) 和式 (2) 写成离差的形式：
将式 (3) 左右两边同时乘以 ，并对 求和可得：
简单整理之后可得：
因此，
其中， 实际上是 对 回归的系数，刻画了 和 之间的相关系数，即：
因此遗漏变量偏误的表达式为：
等式 (4) 的推导利用到 和 的无偏性。此外 和 的方差的表达式为：
其中，，。
这其中有几点值得注意，如果 和 高度相关，但是 对 的影响很小，则 的偏误事实上很小，也就是说 和 之间很接近。此时 。相反，如果 和 之间相关性很小，但是 对 的影响很大，此时 的偏误事实上也很小，但是 。
当然对于上述遗漏变量偏误的推导也有一些别的方法，比如：
不过上述的推导过程尽管看似复杂，但是很容易推广到更多元的情形之下，因此我们在这里将其较为详细地列示出来，当然如果读者掌握更高级的矩阵知识的话，也可以有更为简便的办法，在这里就不再赘述了。
由于在实证分析中，我们很难将全部的解释变量都纳入进模型当中，因此将不可避免地受到遗漏变量的影响。尽管如此，我们仍可以分析结果对于遗漏变量偏误的敏感性。
一些经典的研究遗漏变量对于实证结果的影响的文献 (Altoji，2005；Oster，2019) 均要求外生控制变量假设，即遗漏变量必须与控制变量不相关，而这一假设通常是一个很强而且不太合理的假设。本文主要探讨在 Diegert 等 (2022) 提出的敏感性分析方法，该方法允许内生控制变量的存在。在这一框架下，有两个核心问题：
在放松的假设下， 的一致估计是什么？或者说在备择假设下， 的边界是什么？ 在 的一个假设不被拒绝之前，我们可以在多大程度上放松外生性假设？也就是截断点 (breakdown point)：在一个假设被推翻之前，基准假设的最大可放松程度。我们希望这个程度越大越好，越大表示我们的结论越稳健。
接下来，我们将使用命令 regsensitivity 对上述两种敏感性进行分析。
2. 方法初探
这里采用 Bazzi 等 (2020) 的数据进行演示，我们首先来复现 Bazzi 等 (2020) 中表 3 第 7 列。
. lxhget regsensitivity.pkg, install
. lxhuse bfg2020.dta, clear
. local y avgrep2000to2016
. local x tye_tfe890_500kNI_100_l6
. local w1 log_area_2010 lat lon temp_mean rain_mean elev_mean d_coa d_riv d_lak ave_gyi
. local w0 i.statea
. local w `w1' `w0'
. local SE cluster(km_grid_cel_code)
. reg `y' `x' `w', `SE'
Linear regression Number of obs = 2,036
F(39, 379) = .
Prob > F = .
R-squared = 0.3321
Root MSE = 9.6368
(Std. err. adjusted for 380 clusters in km_grid_cel_code)
------------------------------------------------------------------------------------------
| Robust
avgrep2000to2016 | Coefficient std. err. t P>|t| [95% conf. interval]
-------------------------+----------------------------------------------------------------
tye_tfe890_500kNI_100_l6 | 2.055 0.349 5.88 0.000 1.368 2.741
log_area_2010 | 0.276 0.980 0.28 0.778 -1.651 2.203
lat | 2.265 1.101 2.06 0.040 0.100 4.430
lon | 0.011 0.291 0.04 0.970 -0.562 0.584
temp_mean | 1.627 1.068 1.52 0.128 -0.473 3.728
rain_mean | 0.016 0.005 3.58 0.000 0.007 0.026
elev_mean | 0.015 0.004 4.10 0.000 0.008 0.023
d_coa | 0.000 0.000 2.62 0.009 0.000 0.000
d_riv | 0.000 0.000 3.10 0.002 0.000 0.000
d_lak | 0.000 0.000 0.07 0.945 -0.000 0.000
ave_gyi | -3.780 10.810 -0.35 0.727 -25.035 17.475
statea |
5 | -4.214 3.386 -1.24 0.214 -10.872 2.445
8 | -27.317 6.247 -4.37 0.000 -39.600 -15.034
12 | 4.628 3.355 1.38 0.169 -1.968 11.224
13 | 0.540 2.644 0.20 0.838 -4.658 5.738
17 | -10.258 3.787 -2.71 0.007 -17.705 -2.811
18 | -5.924 3.497 -1.69 0.091 -12.801 0.952
19 | -18.027 4.514 -3.99 0.000 -26.903 -9.151
20 | 1.599 4.634 0.35 0.730 -7.512 10.710
21 | -0.504 3.186 -0.16 0.874 -6.768 5.760
22 | 0.894 3.277 0.27 0.785 -5.549 7.337
26 | -14.063 4.406 -3.19 0.002 -22.725 -5.401
27 | -18.103 4.821 -3.75 0.000 -27.583 -8.623
28 | -6.931 3.676 -1.89 0.060 -14.159 0.297
29 | -4.170 3.902 -1.07 0.286 -11.843 3.502
31 | -1.343 4.738 -0.28 0.777 -10.658 7.972
35 | -40.780 9.264 -4.40 0.000 -58.996 -22.564
36 | -9.822 4.689 -2.09 0.037 -19.041 -0.602
37 | -13.538 4.242 -3.19 0.002 -21.878 -5.197
38 | -11.982 5.512 -2.17 0.030 -22.821 -1.143
39 | -6.191 3.656 -1.69 0.091 -13.378 0.997
40 | 13.600 4.881 2.79 0.006 4.004 23.197
42 | -3.146 4.426 -0.71 0.478 -11.850 5.557
46 | -11.847 5.102 -2.32 0.021 -21.878 -1.816
47 | -3.541 2.794 -1.27 0.206 -9.035 1.953
48 | 12.826 4.174 3.07 0.002 4.619 21.033
51 | -0.812 4.048 -0.20 0.841 -8.771 7.147
54 | -3.244 3.570 -0.91 0.364 -10.264 3.776
55 | -18.929 4.504 -4.20 0.000 -27.785 -10.073
56 | -19.023 9.729 -1.96 0.051 -38.153 0.107
_cons | -73.535 57.847 -1.27 0.204 -187.277 40.206
------------------------------------------------------------------------------------------接着，进行默认的敏感性分析。
. regsensitivity `y' `x' `w', compare(`w1')
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 80.4%
Other Params : cbar = 1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.05, 2.05 ]
0.095 [ 1.91, 2.20 ]
0.196 [ 1.76, 2.35 ]
0.296 [ 1.59, 2.52 ]
0.397 [ 1.41, 2.70 ]
0.497 [ 1.20, 2.91 ]
0.592 [ 0.95, 3.16 ]
0.693 [ 0.61, 3.50 ]
0.793 [ 0.07, 4.04 ]
0.894 [ -1.05, 5.16 ]
0.989 [-58.90, 63.01 ]
0.989 [ -inf, +inf ]
--------------------------------------------------------------------------------事实上，上述结果已经回答了我们在引言中提到的两个问题： 在一系列假设之下的边界在哪里； 在何时会被推翻。下面我们分别对于这两点展开讨论。
3. 边界
考虑如下模型：
其中， 是可观测的， 是遗漏变量，可能与 相关。注意此处 的记号，我们下文还会提到 和 。现将其含义总结于下表之中：
| 符号 | 含义 |
|---|---|
| 对 的回归中 的系数 | |
| 对 回归中 的系数 | |
| 对 回归中 的系数 |
的联合分布约束是由三个标量敏感性系数 所控制。给定可观测变量 的联合分布以及敏感性参数值，Diegert 等 (2022) 介绍了如何计算 在可识别集合中的的上下边界 。
所谓可识别集合是指在现有假设下给定可观测数据的分布， 可以被一致估计的值。当 且 时 无法被点估计，因此我们需要估计其边界。在 Stata 中要计算 可以采用 regsensitivity bounds 的子命令完成。
regsensitivity bounds depvar indepvar controls, options...其中，depvar 是被解释变量 ，indepvar controls 是解释变量 。值得注意的是，在 regsensitivity 中，解释变量的顺序至关重要。第一个解释变量 是我们感兴趣的变量，也是要进行敏感性分析的变量，而控制变量则是一些我们并不关心的变量。
默认情况下，regsensitivity bounds 计算在 和 保持不变的前提下，对于 的一系列值，对应的边界。默认情况下 ，，为指定不同的 可以使用 cbar 的选项，例如：
. regsensitivity bounds `y' `x' `w', compare(`w1') cbar(.1)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 119%
Other Params : cbar = .1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.05, 2.05 ]
0.394 [ 1.45, 2.66 ]
0.808 [ 0.74, 3.37 ]
1.202 [ -0.02, 4.13 ]
1.617 [ -0.93, 5.04 ]
2.032 [ -2.02, 6.13 ]
2.425 [ -3.32, 7.43 ]
2.840 [ -5.17, 9.28 ]
3.234 [ -7.86, 11.97 ]
3.648 [-13.63, 17.74 ]
4.042 [-75.24, 79.35 ]
4.063 [ -inf, +inf ]
--------------------------------------------------------------------------------为将结果绘制出来，可以使用命令 regsensitivity plot。
注意在上述命令中，我们加入了 compare(varlist) 的选项，这一选项是在指明哪些控制变量被放入 而不是 . 这些控制变量被称作比较控制变量 (comparison controls), 他们是用来校准敏感性参数 的。通过将更多的变量加入比较控制变量集中后，在给定敏感性参数时， 的可识别集合通常而言会变大。而当这一选项被省略掉之后，全部的控制变量都将被包括进 中，比如说
. regsensitivity bounds `y' `x' `w', cbar(.1)
Regression Sensitivity Analysis, Bounds
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.708
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.027
R2(medium) = 0.332
Var(Y) = 136.320
Var(X) = 1.257
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 29.7%
Other Params : cbar = .1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.05, 2.05 ]
0.128 [ 1.20, 2.91 ]
0.262 [ 0.25, 3.86 ]
0.396 [ -0.77, 4.88 ]
0.531 [ -1.91, 6.02 ]
0.665 [ -3.24, 7.35 ]
0.800 [ -4.88, 8.99 ]
0.934 [ -7.07, 11.18 ]
1.069 [ -10.45, 14.56 ]
1.203 [ -17.46, 21.57 ]
1.331 [-106.48, 110.59 ]
1.336 [ -inf, +inf ]
--------------------------------------------------------------------------------
. regsensitivity plot不难发现，当全部的控制变量都被包含进 和 时，可识别集合就变成了全体实数集 , 而当 中去除州层面的固定效应之后，只有当 时，可识别集合才会变为全体实数集 。为直接比较两种 选择下的边界，我们可以手动地设置 rxbar 值，使其在两种情形下是相同的。
由于 regsensitivity bounds 的回归结果被存储在 e(idset_table) 中，因此我们可以提取其中的 rxbar 的值，然后将其作为选项添加到从 中去除州层面的控制变量之后的敏感性回归命令中，即：
. forvalues i=1/12{
1. local rxbar `rxbar' `= e(idset_table)[`i', 1]'
2. }
. regsensitivity bounds `y' `x' `w', compare(`w1') cbar(.1) rxbar(`rxbar')
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 119%
Other Params : cbar = .1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.0548, 2.0548 ]
0.128 [ 1.8628, 2.2468 ]
0.262 [ 1.6550, 2.4546 ]
0.396 [ 1.4408, 2.6687 ]
0.531 [ 1.2198, 2.8898 ]
0.665 [ 0.9911, 3.1184 ]
0.800 [ 0.7540, 3.3555 ]
0.934 [ 0.5078, 3.6017 ]
1.069 [ 0.2513, 3.8583 ]
1.203 [-0.0166, 4.1262 ]
1.331 [-0.2829, 4.3924 ]
1.336 [-0.2937, 4.4032 ]
--------------------------------------------------------------------------------不难看出，当全部的控制变量都被包括进 时，每一个相同的 对应的边界要更窄。为了比较 的不同值，我们可以将一个数列放进 cbar 的选项当中。此时 regsensitivity bounds 命令将会绘制一张图，而不是以表格的形式呈现在命令窗口
. regsensitivity bounds `y' `x' `w',compare(`w1') cbar(0(.2)1)在默认情况下，regsensitivity plot 将绘制到对于每一个 可识别集合变为全体实数集 。当然我们也可以限制 的取值范围来考察可识别集合与 0 相交的点。
4. 截断前沿
regsensitivity bounds 的结果显示了对于一个给定假设下，参数 的截断点。对于一个假设 ，截断点是敏感性参数 集合的下界。在 这个集合中， 这一假设对于可识别集合中的每一个 都不成立，正式的定义如下：
regsensitivity 可以处理的假设形式为：对于任意的 ，。默认情况下，假设形式为 。其中 是 对 回归中 的系数。在这个例子中由于 ，因此默认情况下的假设为 。
regsensitivity bounds 显示的是去除掉州层面的固定效应的 ，。为进一步看清截断点是如何随着敏感性参数 的变化而变化的，我们可以使用 breakdown 的子命令。例如：
. regsensitivity breakdown `y' `x' `w', compare(`w1') cbar(0(.1)1)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Hypothesis : Beta > 0 Var(X_Residual) = 0.882
Other Params : rybar = +inf
--------------------------------------------------------------------------------
cbar rxbar(Breakdown)
--------------------------------------------------------------------------------
0.000 135.0%
0.100 119.5%
0.200 108.0%
0.300 99.3 %
0.400 92.7 %
0.500 87.6 %
0.600 83.9 %
0.700 81.4 %
0.800 80.4 %
0.900 80.4 %
1.000 80.4 %
--------------------------------------------------------------------------------
. regsensitivity plot为了检验给定某个 时， 这样的假设，我们可以指明 beta(b,lb)。其中，lb 是 lower bound 的意思。选项 beta 也可以接受一个数列以检验一系列的假设。例如，对于一系列的 值，我们想检验 ，可以输入以下代码：
. regsensitivity breakdown `y' `x' `w', compare(`w1') beta(-1(.2)1 lb)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Hypothesis : Beta > Beta(Hypothesis) Var(X_Residual) = 0.882
Other Params : cbar = 1, rybar = +inf
--------------------------------------------------------------------------------
Beta(Hypothesis) rxbar(Breakdown)
--------------------------------------------------------------------------------
-1.000 89.1 %
-0.800 87.9 %
-0.600 86.5 %
-0.400 84.8 %
-0.200 82.8 %
0.000 80.4 %
0.200 77.4 %
0.400 73.8 %
0.600 69.5 %
0.800 64.1 %
1.000 57.5 %
--------------------------------------------------------------------------------当然，我们也可以通过指定 beta(b ub)来检验 。其中，ub 表示 upper bound。例如，我们想检验对于一系列的 ， ，可以输入以下代码：
. regsensitivity breakdown `y' `x' `w', compare(`w1') cbar(0(.1)1) beta(4 ub)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Hypothesis : Beta < 4 Var(X_Residual) = 0.882
Other Params : rybar = +inf
--------------------------------------------------------------------------------
cbar rxbar(Breakdown)
--------------------------------------------------------------------------------
0.000 128.1%
0.100 114.0%
0.200 103.6%
0.300 95.7 %
0.400 89.6 %
0.500 85.0 %
0.600 81.7 %
0.700 79.5 %
0.800 78.8 %
0.900 78.8 %
1.000 78.8 %
--------------------------------------------------------------------------------为了将上表可视化，我们可以结合 regsensitivity bounds 的命令绘制可识别集合和 4 的交点。
5. 回归结果的解读
无论是 regsensitivity bounds 还是 regsensitivity breakdown，他们的输出结果中都包含着一个描述性统计表，其含义如下：
Number of observations：观测值个数；Beta(short)： 对 回归中 的系数；Beta(medium)： 对 回归中 的系数；R2(short)： 对 回归的 ；R2(short)： 对 回归的 ；Var(Y)：的方差；Var(X)： 的方差；Var(X_Residual): 的方差， 表示 对 回归的残差。
值得注意是，在这张回归表格当中， 均表示 ，其含义同类似。
6. 参考文献
Masten M A, Poirier A. Inference on breakdown frontiers[J]. Quantitative Economics, 2020, 11(1): 41-111. -PDF- Bazzi S, Fiszbein M, Gebresilasse M. Frontier culture: The roots and persistence of “rugged individualism” in the United States[J]. Econometrica, 2020, 88(6): 2329-2368. -PDF- -Replication- Oster E. Unobservable selection and coefficient stability: Theory and evidence[J]. Journal of Business & Economic Statistics, 2019, 37(2): 187-204. -PDF- Altonji J G, Elder T E, Taber C R. Selection on observed and unobserved variables: Assessing the effectiveness of Catholic schools[J]. Journal of political economy, 2005, 113(1): 151-184. -PDF- Altonji J G, Elder T E, Taber C R. Using selection on observed variables to assess bias from unobservables when evaluating swan-ganz catheterization[J]. American Economic Review, 2008, 98(2): 345-50. -PDF- Chalak K. Identification of average effects under magnitude and sign restrictions on confounding[J]. Quantitative Economics, 2019, 10(4): 1619-1657. -PDF- Masten M A, Poirier A, Zhang L. Assessing sensitivity to unconfoundedness: Estimation and inference[J]. arXiv preprint arXiv:2012.15716, 2020. -PDF-
7. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 控制变量, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Stata：控制变量组合的筛选-tuples Stata新命令-pdslasso：众多控制变量和工具变量如何挑选？ 专题：回归分析 调节效应是否需要考虑对控制变量交乘？ 控制变量！控制变量！ 不用太关心控制变量，真的！ 加入控制变量后结果悲催了！ 专题：IV-GMM Lasso一下：再多的控制变量和工具变量我也不怕-T217 专题：断点回归RDD RDD：断点回归可以加入控制变量吗？ Stata：RDD-中可以加入控制变量 专题：其它 锚定情境法（一）：有效控制变量自评偏差
课程推荐：因果推断实用计量方法
主讲老师：丘嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。